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Abstract 

We present a Markov-chain Monte Carlo algorithm of worm type that correctly 
simulates the fully-packed loop model with n = 1 on the honeycomb lattice, and we 
prove that it is ergodic and has uniform stationary distribution. The honeycomb- 
lattice fully-packed loop model with n = 1 is equivalent to the zero-temperature 
triangular-lattice antiferromagnetic Ising model, which is fully frustrated and noto- 
riously difficult to simulate. We test this worm algorithm numerically and estimate 
the dynamic exponent Zexp = 0.515(8). We also measure several static quantities of 
interest, including loop- length and face-size moments. It appears numerically that 
the face-size moments are governed by the magnetic dimension for percolation. 
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1 Introduction 



The antiferromagnetic Ising model on the triangular lattice is of long-standing 
interest since it provides a canonical example of geometric frustration: it is 
topologically impossible to simultaneously minimize the interaction energies 
of all three edges of an elementary triangular face. Recall that the Ising model 
on finite graph G = (V, E) is defined by the measure 

/XG,/5(a)oce-^^('^), cre{-l,+lV, (1.1) 

where the Hamiltonian, H, for the zero-field nearest-neighbor Ising model is 
simply 

Hia) = -Y: ^.^r (1-2) 

ij£E 

The coupling f3 > {(3 < 0) corresponds to a ferromagnetic (antiferromag- 
netic) interaction. We say an edge ij is satisfied if the Ising interaction en- 
ergy of its two endpoints, —Pa^aj, is minimized. An edge is therefore satisfied 
if its endpoints have parallel (anti-parallel) Ising spins in the ferromagnetic 
(antiferromagnetic) case. It is clear that a given elementary face of the tri- 
angular lattice can have at most two satisfied edges in an antiferromagnetic 
Ising model. The zero-field triangular-lattice antiferromagnetic Ising model 
has an exponentially large number of ground state degeneracies, leading to 
non- vanishing entropy per spin pQ. Although the model is disordered at all fi- 
nite temperatures, at zero temperature the two-point correlation function de- 
cays algebraically [2j, and so the model has a zero-temperature critical point. 
By generalizing (11. 2p to include anisotropic couplings, dilution, longer range 
interactions, or a (staggered) magnetic field, rich phase diagrams have been 
observed 



Frustrated systems are notoriously difficult to simulate. Naive algorithms 
such as single-spin-flip dynamics are inefficient at low temperatures, and be- 
come non-ergodi(0 at zero temperature. Indeed, even the best cluster algo- 
rithms [T0|ll|ll2j for simulating low temperature frustrated Ising models be- 
come non-ergodic at zero temperature, although ergodicity can supposedly 
be obtained by augmenting the cluster dynamics with single-spin-flip dynam- 
ics [TOllTl] (a similar hybrid approach is applied to the string dynamics dis- 
cussed in [13]). Cluster algorithms have deflned the dominant paradigm for 
efficient Monte Carlo simulations of critical lattice models ever since the sem- 
inal work of Swendsen and Wang [14J. A more recent idea which is showing 
great promise however is the idea of worm algorithms, first discussed in the 



^ Following the typical usage in the physics literature, we take ergodic as synony- 
mous with irreducible. Recall that a Markov chain is irreducible if for each pair of 
states i and j there is a positive probability that starting in i we eventually visit j, 
and vice versa. 
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context of classical spin models in [T3] (see also [IS])- The key idea behind 
the worm algorithm is to simulate the high-temperature graphs of the spin 
model, considered as a statistical- mechanical model in their own right. The 
worm algorithm for the Ising model was recently studied in some detail in [17], 
and it was observed to possess some unusual dynamic features (see also [18]). 
Indeed, despite its local nature, the worm algorithm was shown to be extraor- 
dinarily efficient - comparable to or better than the Swendsen-Wang (SW) 
method - for simulating some aspects of the critical three-dimensional Ising 
modelFI. Given this success, a natural question to ask is whether one can de- 
vise a valid worm algorithm to simulate a fully frustrated model such as the 
zero-temperature triangular-lattice antiferromagnetic Ising model. The short 
answer is yes, as we demonstrate in Section [21 although it does require a little 
thought. 

Worm algorithms provide a natural way to simulate Eulerian- subgraph models. 
Given a finite graph G = {V, E), we call a bond configuration A (1 E Eulerian 
if every vertex in the subgraph {V,A) has even degree (i.e., every vertex has 
an even number of incident bonds; zero is allowed). The set of all such bond 
configurations defines the cycle space of G, denoted C{G). Perhaps the simplest 
class of Eulerian-subgraph model is defined on the cycle space of a finite graph 
G = (V, E), for n,w > 0, by the probability measure 

0G,«,,n(^) oc n<^^w\^\, A e C(G), (1.3) 

where c{A) is the cyclomatic number of the spanning subgraph {V,A). Note 
that on graphs of maximum degree < 3 the only possible Eulerian subgraphs 
consist of a collection of disjoint cycles, or loops, and c{A) is then simply the 
number of such loops. Consequently, Eulerian-subgraph models often go by 
the name of loop models, and the honeycomb lattice, being a 3-regular graph, 
has played a distinguished role in the literature on such loop models P2 23]. 
These geometric models play a major role in recent developments of conformal 
field theory [21] via their connection with Schramm Loewner evolution (SLE) 



In this work we focus on the case n = 1, and we write (pCw '■= 4>g.w,i- In this 
case it can be seen that (11.31) corresponds to an Ising model on G, and also 
to an Ising model on the dual graph G*, when G is planar. Indeed, it is an 

^ In our opinion, the conventional wisdom that local algorithms are a priori less 
efficient than cluster algorithms does not bear scrutiny. Both the worm algorithm 
and the Sweeny algorithm |19j (a local algorithm for simulating the random-cluster 
model) have efficiencies |17|20) comparable to, and in some instances better than, 
cluster algorithms such as SW or the Chayes-Machta algorithm [21]. Indeed, both 
algorithms can display critical speeding-up [20] in certain situations, (although ad- 
mittedly the Sweeny algorithm suffers from some algorithmic complications due to 
the need for efficient cluster-finding subroutines). 
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elementary exercise to derive the following two identities relating the partition 
functions of the Ising and Eulerian-subgraph models 



'Ising 

G,/3 — 



(2l^l coshl^l /?) Z, 

/<-) yEulerian 




(1.4) 
(1.5) 



'Ising 

G*,f3 — 



Note that (11.41) corresponds to an Eulerian-subgraph model with positive 
weights only when /3 > 0, i.e. in the ferromagnetic case, and that only the re- 
gion < w < 1 is covered by the correspondence. By contrast, (11.51) gives pos- 
itive weights for all /5 G M and corresponds to the whole region < w < +oo. 
Since tanh(/5) is small when (3 is small the relation (11.41) is commonly referred 
to as the high-temperature expansion of the Ising model. Similarly, since e~'^^ 
is small when (5 is large, i.e. in the ferromagnetic regime at low temperatures, 
(II. 5p is commonly referred to as the low temperature expansion of the Ising 
model. However since e~'^^ is in fact large in the antiferromagnetic regime at 
low temperatures we shall refrain from using this terminology. The n = 1 
Eulerian-subgraph model with < w < 1 thus corresponds to ferromagnetic 
Ising models on both G and G*, while the model with 1 < w < +oo cor- 
responds to an antiferromagnetic Ising model on G* . In particular, we point 
out that the honeycomb-lattice loop model with w > 1 corresponds to the 
triangular-lattice antiferromagnetic Ising model. 

The honeycomb-lattice loop model exhibits an interesting phase diagram that 
has been studied in detail [2511^ . When n = 1 it undergoes an Ising phase 
transition at Wc = l/^/S from a disordered phase when w < Wc to a densely- 
packed phase when Wc < w < +oo. Interestingly, the entire region Wc < w < 
+00 displays critical behavior, and the model is in the two-dimensional per- 
colation universality class. The case w = +oo is of especial interest, and is the 
focus of this article. In this case (pCw is simply uniform measure on the set of 
fully-packed configmations, i.e. the set of all Eulerian subgraphs with the maxi- 
mum possible number of edges, and the model is referred to as the fully-packed 
loop (FPL) model. It should be emphasized that the FPL model is critical, 
although it is in a different universality class to the densely-packed phase [29] . 
Every fully-packed bond configuration is such that each vertex is visited by 
precisely one loop, i.e. each vertex has degree 2. Therefore only two thirds of 
the edges are occupied in any given fully-packed configuration; this fact can 
be seen as a symptom of the frustration of the triangular-lattice antiferro- 
magnetic Ising model. Indeed, according to (11.51) the FPL model corresponds 
to the triangular-lattice antiferromagnetic Ising model at zero temperature. 
Fig. [T] shows a typical fully-packed configuration. 

The essence of the worm idea is to enlarge a configuration space of Eulerian 
bond configurations to include a pair of defects (i.e., vertices of odd degree), 
and then to move these defects via random walk. When the two defects co- 
incide, the configuration becomes Eulerian once more. In the standard worm 
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Fig. 1. Typical fully-packed configuration on the honeycomb lattice with periodic 
boundary conditions. Thick lines denote occupied edges, thin lines denote vacant 
edges. 
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algorithm we view the simulation as a simulation of high-temperature graphs 
of the Ising model on G defined by (11. 4p . This interpretation is only valid for 
Q < w < 1. However, another useful interpretation when G is planar is that 
the worm algorithm simulates an Ising model on G*, and this interpretation is 
valid for all > 0. We shall return to this point in some detail in Section ^I7I\ 
We wish to emphasize here however that if we have a worm algorithm to 
simulate the FPL model, we immediately have an algorithm to simulate the 
zero-temperature triangular-lattice antiferromagnetic Ising model. 

Unfortunately, devising a valid worm algorithm for simulating the FPL model 
is not simply a case of taking the large w limit of the "standard" version of the 
worm algorithm, as presented in [17] . Indeed, as the bond weight w increases, 
the efficiency of the worm algorithm presented in Refs. P^IITT] drops rapidly, 
because the random walker moves ever more slowly. In the limit w ^ oo the 
random walk becomes completely frozen, and the standard worm algorithm 
becomes invalid (the details will be explained in Section [2]). In this work, we 
present a variation of the worm algorithm presented in [l7j which efficiently 
simulates the honeycomb-lattice FPL model, when n = 1. Importantly, we 
prove rigorously that this algorithm is ergodic, and has uniform stationary 
distribution, on the fully-packed configurations. We have tested this worm 
algorithm numerically, and we estimate the dynamic exponent Zexp = 0.515(8). 
(See Section [33] for a precise definition of Zexp-) 

The organization of the current work is as follows. Section [2] reviews the stan- 
dard worm algorithm [T5][T7] for Ising high-temperature graphs, and then in- 
troduces a version to simulate the FPL model. In Section [3 we present the 
results of our simulations of the FPL model using the worm algorithm dis- 
cussed in Section |2] Finally, Section H] contains a discussion. 
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2 Worm algorithms 



We begin witli a review of the standard worm algorithm defined on an arbitrary 
graph, which essentially follows the presentation in and then go on to 
discuss its relationship to the Eulerian-subgraph model on G as well as Ising 
models on G and G*. After demonstrating why the standard version becomes 
non-ergodic as w —>■ oOjWe then present a valid worm algorithm for simulating 
the honeycomb-lattice FPL model. 

2.1 The "standard" worm algorithm 

Fix a finite graph G = {V,E), and for any A C E let dA C V denote the set of 
all vertices which have odd degree in the spanning subgraph {V,A). Loosely, 
dA is just the set of sites that touch an odd number of the bonds in the bond 
configuration A. If m, w e V are distinct we write 

Su,v:={ACE : dA = {u,v}}, 

and 

Sy,v ■.= {ACE : dA = 0}. 

We emphasize that Sy^y = C{G) for every v E V . We take the state space of 
the worm algorithm to be 

S := {(A, M,t>) : -u, f G y and A G Su,v}, 

i.e., all ordered triples (A, m, v) with A <Z E and u,v E V, such that A G Su,v 
Note that if {A, u,v) E S then A is Eulerian iS u = v. Thus the bond configu- 
rations allowed in the state space of the worm algorithm constitute a superset 
of the Eulerian configurations. Finally, we assign probabilities to the configu- 
rations in S according to 

^^{AjUjV) (X dudyW^'^^ {A,u,v)eS, (2.1) 

where dy denotes the degree in G of f G V^. In the following, when we wish 
to refer to the degree of f G K in the spanning subgraph {V, A) we will write 
dy{A). Loosely, dy{A) is simply the number of bonds that touch v in the bond 
configuration A. In this notation we have dy = dy{E). 

The first step in constructing the standard worm algorithm is to consider the 
worm proposal matrix, P^^\ which is defined for all uu' G E and v E V hy 

p(°)p,M,t;) {AAuu',u',v)] = P^''\iA,v,u) iAAuu',v,u')] = -i-, 

(2.2) 
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all other entries being zero. Here A denotes symmetric difference, i.e. delete 
the bond uu' from A if it is present, or insert it if it is absent. It is easy to see 
that P(0) is an irreducible transition matrix on S. According to (12.21) the moves 
proposed by the worm algorithm are as follows: Pick uniformly at random one 
of the two defects (say, v) and one of the edges emanating from v (say, ww'), 
then move from the current configuration (A, m, v) to the new configuration 
{A A vv',u,v'). 



Now we simply use the usual Metropolis-Hastings prescription (see e.g. 
§4]) to assign acceptance probabilities to the moves proposed by P^^\ so that 
the resulting transition matrix, P^, is in detailed balance with (12. ip . Explicitly, 
for all uu' G E and v & V we have 

P^[{A,u,v) ^ {AAuu',u',v)] = Py,[{A,v,u) {AAuu',v,u')] 

1 iF{w) uu'^A (2.3) 
~2^|F(1/w) uu'eA 



where F : [0, +oo] [0, 1] is any function satisfying 

F{z) = zF{l/z) for all z. 



(2.4) 



Two concrete examples of such F which are commonly used in practice are 
F{z) = min(l, z) and F{z) = zj (I + 2;). For a given choice of F, the transitions 
(12. 3p define P^ uniquely since all other transitions occur with zero probability 
except the identity transitions (A, M,f) — > (A, M,t>), whose transition proba- 
bilities are fixed by normalization to be 



Pj,{A,u,v) ^ {A,u,v)\ = \-F{w) 



1 - 



2d„ 



2d, 



(2.5) 



2du 



2d^ 



For any choice of F, one can easily verify that and tt^ are in detailed 
balance. 



2.2 Relation to Eulerian- subgraph and Ising models 



A natural question to ask at this stage is what precisely is the relationship 
between (pcw and the worm transition matrix (12.31) ? To address this question, 
let us consider the Markov chain induced on the subset 

S := {{A,v,v) eS} CS, (2.6) 

in which the bond configurations are Eulerian. More precisely, let's suppose 
that we only observe the worm chain when it is in a state in 5*. This de- 
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fines a new Markov chain, a single step of wfiicli corresponds in tlie old chain 
to the transition (not necessarily in one step) from a state {A, v, v) to an- 
other state {A\v',v'). The new transition probability to move from {A,v,v) 
to {A\v',v') is found by computing the probability that the original chain 
starting in (A, f,w) hits 5* for the first time at state {A\v',v'). This is the 
probability that the chain goes from (A, f , v) to (A', f ', v') in one step (which 
is zero unless A = A' and v = v'), plus the probability that it goes to a state 
outside S and then re-enters S for the first time at (A', t>', v'). A nice discussion 
of this general problem can be found in [31l §6.1], including a proof of 

Lemma 2.1. Let P he an irreducible transition matrix on a finite state space 
S with stationary distribution tt. Define a new Markov chain by only observing 
the original chain corresponding to P when it visits a state in S <Z S. The new 
chain is an irreducible Markov chain on S with stationary distribution 

-Ks = — , s e S. 



As a consequence of Lemma 12.11 the worm Markov chain restricted to the 
Eulerian subspace (12.61) has a stationary distribution given explicitly by 

7rUAv,v) = ( ] 0G,.(A). (2.7) 

Consequently we have 

(X)^^ = (X)<^^,„ (2.8) 

for any observable X : C{G) — ^> M of the original Eulerian-subgraph model, 
and hence we can indeed use the worm algorithm to simulate 4>g,w 

We note that when G is planar (12.81) also implies that the worm algorithm 
correctly simulates the Ising model on G* considered in (11.51) . Indeed, sup- 
pose that G is planar with dual G* = {V*,E*), and consider the two-to-one 
correspondence a Acr from { — 1,+1}^* C{G) where 

A„:={ij eE ■ ai, ^aj,}. (2.9) 

In words, for any spin configuration on G* we draw on G the boundaries of the 
spin domains. It is an elementary exercise to show that for all a G { — 1, +1}^* 
we have 

0G,«,(^<x) = 2^G*,/3(^), w = e-^^, (2.10) 
where /xg*,/? is the mass function of the Ising model on G*, as defined in (11. ip . 
We emphasize that although (I2.10p is often called a low temperature represen- 
tation, it is an exact result valid for all — oo < j3 < +oo, or equivalently for 
all < w < +00. From (12.101) we see explicitly that for any Ising observable 
Y : { — 1, +1}^* M that is even under global spin flips (which is the case for 
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all observables of physical interest in zero field) we have (Y)^^^ = (X)^^^ 
where w = e-^^ and X : C(G) ^ R is defined by X{A„) = Y(a) = Y{-a). 
Consequently f l2.8p does indeed allow us to simulate the Ising model on G* 
using the worm algorithm. 

We should also mention that when w < 1 the worm algorithm can be used to 
simulate properties related to the two-point correlation function of the Ising 
model on G defined by (11.41) . Indeed, it is straightforward to generalize fll.4l) 
to obtain an expansion for the two-point correlation function 

4T(^n^^)Mc,. = E (tanh/3)l^l. (2.11) 

,v 

As an example of the use of (12.111) . consider the observable T>o on S defined 
so that 

VoiA,u,v)=5u,,. (2.12) 
In other words T>q is the indicator for being in 5*. It is straightforward to show 
that provided G is regular we have 

where = J2v&v '^v is the magnetization and the vr^ expectations use w = 
tanh p. In particular, in a translationally invariant system 

= 1/Xg7- (2-13) 

Thus when w < 1 the worm algorithm simulates both an Ising model on G 
and an Ising model on G*. Quantities like Pq depend on the full Markov chain 
on S, and so if one's interest is to obtain quantities related to the two-point 
function for the Ising model on G with w = tanh(/5) then one must consider 
the full Markov chain. However, if one's interest is to compute properties of 
the Eulerian-subgraph model (11.31) . or the corresponding Ising model on G* 
with P = e~^^, then one is only interested in the Markov chain induced on 
S. It is the latter models that are our interest in the present work, and we 
emphasize that in this case the restriction w < 1 does not apply. 

We note, finally, that [32j uses ideas similar to Lemma 12.11 and (12.101) to 
simulate a low temperature Ising spin glass with a worm algorithm. 

2.3 Periodic boundary conditions 

For completeness, we now briefly address the question of the effect of boundary 
conditions when G is a regular lattice. To illustrate, we consider G = M., where 
H denotes a finite subgraph of the honeycomb lattice drawn on a torus as in 
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Fig. [TJ The periodic boundary conditions imply that EI is non-planar, however 
we can still construct the dual lattice T in the usual way, and it is easy to see 
that T is simply a finite piece of the triangular lattice also drawn on a torus. 
It is now no longer the case however that every A G C{M.) defines the domain 
boundaries of an Ising spin configuration; indeed Fig. [1] provides an example 
for which no consistent assignment of Ising spins is possible. Suppose however 
that we let C^(HI) denote the set of all A G C{M) which wind the torus an 
even number of times in both directions. For these configurations there is no 
ambiguity in assigning Ising configurations according to the correspondence 
(12. 9p . and it is easy to see that (12.91) defines a two-to-one correspondence from 
{— 1, +1}^'^^) onto C"''(EI). It is easy to generalize (l2.1Up to show that it is now 
replaced by 

In addition, if one applies Lemma [2. II to the subspace C^(H) then we obtain 

1 tyl^l 

7rMAAv,v) = -- ^, for all A G C+(H). (2.15) 

y z^a'gc+(h) I 

Combining (I2.14p and (I2.15P we see immediately that Wj;i^^-2/3{Aa,v,v) = 
(2/y)/ix,/3(o"). Therefore if we simulate a worm chain on EI with coupling e~^^ 
and only measure this chain when it is both Eulerian and winds the torus an 
even number of times, then we are effectively simulating the Ising model on 
T at inverse temperature /3. 



2.4 A worm algorithm for the honeycomb lattice FPL model 



Thus far we have glossed over an important issue, namely the irreducibility 
of the worm transition matrix P^. It is not hard to see that is irreducible 
whenever F{w) and F{l/w) are both strictly positive. Problems arise as w —>■ 
00 however, since it is easy to show that if F : [0, +00] [0, 1] satisfies (12.41) 
then -F(O) = 0. Consequently, as w —>■ 00 the probabilities for transitions 
that remove an edge vanish. Indeed, all states [A, u,v) E S for which both 
du{A) = du and dv{A) = d^ become absorbing as w ^ 00. This is easy to see 
from (12. 5p since such states have Pyo[{A,u,v) {A,u,v)] = 1 — F{l/w). 

Suppose now that G is fc-regular, i.e. all vertices have degree k. Then {A, u, v) 
will be absorbing when w = +00 iff du{A) = dv{A) = k. Recall that if 
{A,u,v) G S then when u = v the vertex degree du{A) = d^^A) is even, 
whereas when u ^ v both du{A) and d^{A) are odd. Thus if k is even then 
{A, u, v) can be absorbing only if -u = f whereas if k is odd {A, u, v) can be 
absorbing only ii u ^ v. Therefore when k is odd all states {A,v,v) with 
Eulerian A remain non-absorbing; [A, v, v) — >■ [A, v, v) occurs with probabil- 
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ity dy{A)/k < 1 when w = +00. In particular, on the honeycomb lattice we 
can now see that as w ^ 00 all states {A,v,v) with Eulerian A remain non- 
absorbing while all states {A, u, v) with u ^ v and du{A) = dv{A) = 3 become 
absorbing. Therefore once both defects have degree 3 the chain remains in 
that state for eternity. 

How do we resolve this problem? A simple answer is to avoid this trap of end- 
less identity transitions by explicitly forbidding [A, u, v) — > {A, u, v) whenever 
u ^ V. Since, when simulating Eulerian-subgraph models, we only observe the 
chain when it visits an Eulerian state {A, v, v) we may hope that by only mod- 
ifying the transitions from non-Eulerian states we may recover irreducibility 
without sacrificing the correctness of the stationary distribution. We shall see 
that this is indeed possible. 

To this end we now define a new transition matrix, P^, which defines a valid 
Monte Carlo algorithm to simulate the FPL model on the honeycomb lattice, 
i.e. when G = M with HI as defined in Section 12.31 To define the transition 
probabilities P^[{A,v,v) ^ • ] we simply take the limits of (12. 3p 



PL[iAv,v) ^ iAUvv',v\v)] = PL[iA,v,v) iAUvv\v,v')] = -, (2.16) 
and 



All other transitions from {A, v, v) are assigned zero probability; in particular, 
one cannot remove an edge from an Eulerian state. 

To define the transition probabilities P^[{A,u,v) ■] with u v we use 
the following simple rules: first, choose uniformly at random one of the two 
defects, say u. Since u ^ v we must have du{A) G {1,3}. If du{A) = 3 we 
choose uniformly at random one of the three occupied edges incident to u, say 
uu', and we delete it by making the transition {A, u, v) —>■ {A\ uu', u', v). This 
ensures that we can never get stuck when the defects are full - i.e. it removes 
the problem of absorbing states suffered by the w ^ 00 limit of P^. If, on 
the other hand, du{A) = 1 we choose uniformly at random one of the two 
vacant edges incident to u, say uu', and occupy it by making the transition 
{A, u, v) {AUuu', u', v). This guarantees that we cannot produce an isolated 
vertex by moving a degree 1 defect, which is obviously a desirable property 
when one wants to simulate a fully-packed model. These rules correspond to 
the following transition probabilities when u ^ v 



1 



P:^[{A,v,v)-^{A,v,v)] 



dM) 
3 



(2.17) 




(2.18) 
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All other transitions from {A, u, v) with u ^ v are assigned zero probability; 
in particular, no identity transitions are allowed. 

While we hope that the above discussion convinces the reader that provides 
a plausible (and natural) candidate for simulating the FPL model on the 
honeycomb lattice, we of course do not claim that it proves such an assertion. 
A proof of the validity of is presented in Section 12.51 

In terms of a Monte Carlo algorithm, P^ corresponds to Algorithm [H The 
abbreviation UAR simply means uniformly at random. 

Algorithm 1 (Honeycomb-lattice fully-packed loop model). 
loop 

Current state is {A, u, v) 
ifu = v then 

Choose, UAR, one of the 3 neighbors of u (say u') 
if uu' ^ A then 

Perform, UAR, either {A, u, u) {A U uu', u', u) or {A, u, u) 
{A U uu', u, u') 
else if uu' E A then 
{A, u, u) {A, u, u) 
end if 
else ifu^v then 

Choose, UAR, one of the 2 defects (say u) 
if du{A) = 3 then 

Choose, UAR, one of the 3 neighbors of u (say u' ) 
{A, u, v) ^ {A \ uu', u', v) 
else if du{A) = 1 then 

Choose, UAR, one of the 2 vacant edges incident to u (say uu') 
{A, u, v) {AU uu', u', v) 
end if 
end if 
end loop 



2.5 Proof of validity of AlgorithmUl 

This Section provides a rigorous proof of the validity of Algorithm [H Readers 
uninterested in such details may simply choose to trust us and skip to the 
next Section. 

Proving validity of Algorithm [T] boils down to showing that P^ is irreducible 
(in a suitable sense) and that it has the right stationary distribution (in a 
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suitable sense). With regard to the latter question we note that 0h,oo(^) = 
lFa{A)/\Fu\, where 

Fe := {A G C(H) : d^{A) = 2 for all v G ^(H)} 

and Ip^ is its indicator. That is, (/)e,oo is just uniform measure on the set Fe 
of fully-packed configurations on H. 

Let us pause to recall some basic background regarding finite Markov chains 
(see e.g. [33IIM] ). Consider then a Markov chain on a finite state space with 
transition matrix P. We say state i communicates with state j, and write 
i — > j, if the chain may ever visit state j with positive probability, having 
started in state i. We say states i and j intercommunicate, and write i ^ j, if 
i ^ j and j — >■ i. A set of states C is called irreducible ii i ^ j for all i,j G C, 
and it is called closed if Pij = for alH G C and j ^ C. A state i is recurrent 
if, with probability 1, the chain eventually returns to i, having started in i; 
and it is transient otherwise. If every state in C is recurrent (transient) we say 
C itself is recurrent (transient). It can be shown that C is recurrent iff it is 
closed. 

Now let us define the subset of states 

n= {{A,u,v) eS : d^{A) =^ for all x} (2.19) 

The set TZ thus consists of all those states with no isolated vertices, and is 
where all the action takes place when considering P^. We emphasize that the 
set of all bond configurations A for which [A, v,v) eTZ corresponds precisely 
with Fh. 

Proposition 2.2. TZ is closed. 

Proposition 2.3. TZ is irreducible and S\TZ is transient. 

Thus when running Algorithm [T] we are free to begin in any state in S, and (due 
to the transience of iS \ 7^) with probability 1 the chain will end up inside 7Z, 
from where (due to TZ being closed) the chain then never leaves. Furthermore, 
(due to the irreversibility of TZ) all states in TZ will eventually be visited. 
Finally, we have the following explicit form for the stationary distribution of 
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Proposition 2.4. The unique stationary distribution of is vr^ where 
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u 7^ 17, + dy{A) = 


4 


>/2 
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7^^, 


U7^v,du{A) + d^,{A) = 
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anc? A zs finite and positive. 

In particular, tt^ is constant on tlie states {A, v, v) G 71. It follows that if 
we consider the Markov chain constructed by only measuring the chain 
when the defects coincide, then Lemma 12.11 implies that P'oo has stationary 
distribution 



n'oo{A,v,v) 



0Fe,oo(^) 



as desired. Consequently {X):^ 
of the FPL model. 



V 



(2.20) 
(2.21) 



for any observable X : Fa ^ '. 
We now conclude this Section with proofs of Propositions 12. 2^ 12. 3^ and 12.41 



Proof of Proposition All the states in 5 \ 7^ have at least one isolated 
vertex, while the states in TZ have none. Since P^ only allows transitions 
that add/remove at most one edge, the only states {A,u,v) G 71 which could 
possibly make a transition to a state with an isolated vertex are those with at 
least one vertex u with du{A) = 1, and the transition would need to remove 
the edge uu' G A. However, we have 

P!^[{A,u,v) ^ {A\uu\u\v)]=0. 

In fact, if du{A) = 1 the only non-zero P^[{A,u,v) ■] that correspond to 
the removal of an edge are of the form 

P^[{A, u, v)^{A\ vv', u, v')] = 1/6 > 

with d^{A) = 3. Therefore the only possible way a transition could remove 
uu' was if t> = m' and we made the transition {A,u,v) {A\ uv,u,u). See 
Fig. O Such a transition would indeed occur with positive probability. However 
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Fig. 2. Here thick lines denote occupied edges, thin hnes denote vacant edges, while 
dashed lines denote edges whose occupation status is undecided. Periodic boundary 
conditions are imposed. A transition capable of creating an isolated vertex can only 
occur from a state for which the neighborhoods of the defects are as shown. There 
clearly exist bond configurations in S with defect neighborhoods as shown, however 
Lemma lA. II implies that no such bond configurations exist in TZ. 
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Lemma lA.ll guarantees that there do not exist any states {A, u,v) & TZ with 
M ~ w and du{A) ^ Therefore 

P:^[{A,u,v)^{A\u',v')] = 

whenever {A, u,v) eTZ and {A, u,v) E S \ 7l. □ 

Proof of Proposition \2.3[ Let H E Fm denote the fully-packed configuration 
in which every horizontal edge is occupied and every vertical edge is vacant. 
We begin by proving that every state in S communicates with {H, w,w) eTZ 
for some w. We make frequent use of the lemmas listed in Appendix [Bi 

Suppose then that {A, u, v) G S. We can generate a new state from {A, u, v) 
via the map f : S ^ S with f{A,u,v) defined by the following prescrip- 
tion: 

if du{A) = 0, 1 then 

Choose a vacant horizontal edge uu' 

return (A U uu', u', v) 
else if du{A) = 3 then 

Choose the occupied vertical edge uu' 

return {A \ uu', u', v) 
else if du{A) = 2 then 

if there are no vacant horizontal edges then 
return {H, u, u) 

1 1 Note that it must be the case that {A, u, v) = {H, u, u) 
else 

Choose a vacant horizontal edge ww' for which {A, u, u) — > {A, w, w) 
/ 1 Lemma [B.2I guarantees that such a ww' exists 
return (A U ww' , w' , w) 
end if 
end if 
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Fig. 3. Example of repeated application of / to a configuration (A,u,v). Starting 
from the initial configuration in (a), application of / removes an occupied vertical 
edge resulting in the configuration f{A,u,v) shown in (b). Continuing in this way, 
alternately adding vacant horizontal edges and removing occupied vertical edges, 
finally results in the configuration f^^{A,u,v) = {H,v,v) shown in (c). Thick lines 
denote occupied edges, thin lines denote vacant edges. Periodic boundary conditions 
are imposed. 

I I I f I I I I 1 I I I I I i I I I 

I I I I I 1 I I I I 



III III I I I 

(a) iA,u,v) (b) f{A,u,v) (c) f\A,u,v) = {H,v,v) 

The first observation to make is that for any {A, u,v) & S we have {A, u, v) — > 
f{A,u,v). Indeed, if du{A) = 0, 1 or 3 we simply have 

P:^[{A,u,v)-^ fiA,u,v)]>0. 

If du{A) = 2 and there are no vacant horizontal edges then it must be the 
case that {A,u,v) = {H,u,u) = f{H,u,u), so {A,u,v) ^ f{A,u,v) follows 
trivially. Finally, if du{A) = 2 and there exists at least one vacant horizontal 
edge then Lemma IB. 21 guarantees that at least one such edge ww' satisfies 
{A,u,u) — >• {A,w,w) and since 

[{A, w,w) ^ {AU ww', w', w)] = 1 /6 

it follows that {A,u,u) — > (A U ww',w',w). So we indeed have {A,u,v) 
f{A,u,v) for any {A,u,v) G S, and in fact {A,u,v) — > f"'{A,u,v) for any 
n eN, where denotes n-fold composition of / with itself, i.e. 

r = /o/o...o/ 

n times. 

Now, whenever {A, u, v) ^ {H, w, w) for some w, the state f{A, u, v) has either 
one less occupied vertical edge, or one more occupied horizontal edge, than 
{A,u,v). Therefore, since there are only a finite number of horizontal and 
vertical edges, if we start in any {A, u,v) E S and apply / repeatedly then we 
must eventually have f"'{A,u,v) = {H,w,w) for some w, with n necessarily 
finite. See for example Fig. El It then immediately follows that {A, u, v) — 
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Suppose now that {A, u,v) E S \ 7l. As we have just shown, there is at least 
one state {H, w,w) E IZ with which (A, m, v) communicates, i.e. (A, u, v) — > 
(if, w, w), and there is thus a non-zero probabihty that starting in (A, u, f ) a 
finite number of transitions will take us to (if, w, w). But since (if, w,w) eTZ 
and Proposition 12.21 tells us that 7^ is closed, there is zero probability of ever 
leaving TZ again, and in particular there is zero probability of ever returning 
to {A, u,v) E S\TZ. There is therefore a non-zero probability that starting in 
{A,u,v) E S \ 7l we never return to {A,u,v). Therefore the state {A,u,v) is 
transient and it follows at once that in fact the whole space iS\7^ is transient. 

Now let us turn our attention to the irreducibility of TZ. It is clear that 
f{A,u,v) E IZ whenever {A,u,v) E 7Z. Furthermore, whenever {A,u,v) E 7Z 
we have f{A, u, v) ^ {A, u, v). To see this we note: we can never have du{A) = 
when {A,u,v) E TZ; if du{A) = 1 then Lemma lB.51 implies {A,u,v) ^ 
f{A,u,v); if du{A) = 3 then Lemma lB.61 implies {A,u,v) ^ f{A,u,v); if 
du{A) = 2 then Lemma lB.31 implies {A,u,u) ^ {A,w,w) for all w, and if 
ww' is vacant Lemma [B.4I implies that {A,w,w) ^ {AU ww',w',w), so that 
{A, u, u) ^ {AU ww', w', w). 

Therefore we now see that for any {A, u,v) eTZ we have {A, u, v) ^ f{A, u, v), 
and indeed {A,u,v) ^ f"^^{A,u,v) for any n eN. Since, as argued above, we 
must have (if, w, w) = f"{A, u, v) for some w and finite n, it immediately fol- 
lows that {A,u,v) {H,w,w). Since every {A,u,v) E TZ intercommunicates 
with (if, w,w) eTZ for some (in fact all) w it follows that TZ is irreducible. □ 

Remark 2.1. The careful reader will notice that there is some ambiguity in the 
definition of / presented in the proof of Proposition 12.31 For instance, if there 
is more than one vacant horizontal edge which one should we choose? Such 
careful readers can easily construct an appropriate rule to make the choice of 
this edge precise (or make the choice of edge random and view / as a random 
variable). The validity of the proof is independent of any such technical details 
and so we have deliberately swept such issues under the proverbial rug. 



Proof of Proposition 2.4' It is straightforward (if a little tedious) to prove that 



is a stationary distribution for P^ by simply considering each of the eight 



cases in the definition of vr^, explicitly computing the right-hand side of 



XA, u,v) = J2 '^'oo{B, X, y)P^[{B, x, y) {A, u, v)] 

(B,x,y)&S 



and verifying that it equals the left-hand side, for every {A, u, v) E S. We omit 
the details. 



Clearly, the constant A appearing in the definition of vr^ must be chosen so 
that Y.{A,u,v)£S — 1; but its exact value is not really of any concern 

to us. We simply observe that it is some well defined finite positive number. 
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Indeed it is elementary to derive the upper and lower bounds 1/A > V^|-Fh| > 
and 1/A < |7^|. 

Since S has only one closed irreducible set of states, 71, it can have only one 
stationary distribution, so tt^ is unique. □ 



3 Numerical results 

We simulated the FPL model on an L x L honeycomb lattice with periodic 
boundary conditions using Algorithm [H We studied fourteen different system 
sizes in the range 6 < L < 900, each being a multiple of 3. 

3.1 Observables measured 

We measured the following observables in our simulations. All observables were 
measured only when the defects coincided, except for Vq which was measured 
every step. 

• The number of loops A/; (cyclomatic number) 

• The mean-square loop length 

£2 := 5] (length of i"" hopf (3.1) 

i=l 

• The sum of the nth powers of the face sizes 

^n:=El/r (3.2) 

/ 

Every A E can be decomposed into a number of faces, each consisting 
of a collection of elementary hexagons, such that every pair of neighboring 
elementary hexagons in H which share an unoccupied edge in A belong to 
the same face. The size |/| of face / is then simply the number of elementary 
hexagons which it contains. We considered n = 2 and n = 4. 

• Vo as defined in fl^IT^ 

From these observables we computed the following quantities: 

• The loop-number density rii := {Ni) / L"^ 

• The loop-number fluctuation Ci := Yai{N'i) / L?' 

• The (normalized) expectation of C2 

L2 ■■= {C,)/L' 
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• The (normalized) expectation of Q2 and Q4 

Y (3.3) 
G, := j^iG,) 

• The ratio Qg := G^Gi 

• The mean number of iterations of the full worm chain between visits to the 
Eulerian subspace 

Te := l/{Vo).'^ 

Remark 3.1. In the case of the FPL model the number of bonds J^{A) = \A\ 
is constant since every vertex has degree 2, and so A/" is a trivial observable in 
this case, unlike the case for Ising high-temperature graphs |17] . 

3.2 Static data 

For each observable O = T^, L2, G2 we performed a least-squares fit to the 
simple finite-size scaling ansatz 

0{L) = L'^-^^°{Oo + CiL^^o-'^ + 021^' + Osiy^). 

The Oi term arises from the regular part of the free energy, while the O2 
and terms correspond to corrections to scaling. The correction-to-scaling 
exponents were fixed to yi = —2 and y2 = —3, and of course c? = 2. As a 
precaution against corrections to scaling we impose a lower cutoff L > L^i^ 
on the data points admitted to the fit, and we studied systematically the 
effects on the fit of varying the value of Lj^in- We estimate 

= 0.2499(2), 
Xl^ = 0.2498(4), (3.4) 
= 0.1040(3). 

According to [29] the magnetic scaling dimension of the n = 1 FPL model is 
Xh = 1/4. From 03.41] we therefore conjecture that in fact 

Xt^ = Xl, =Xn = 1/4. (3.5) 

In particular, we expect the number of iterations of Algorithm [1] between visits 
to the Eulerian subspace to scale like LF'~'^^h = L^/"^, 

We remark that (13. 4p suggests is very close (perhaps equal) to X^'^'^'^ = 
5/48, the magnetic scaling dimension for models in the two-dimensional per- 
colation universality class. Here is a hand-waving argument suggesting that 
in fact Xg2 = ^h^'^'^ might be an identity: For general n it is known [29j that 
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a 




Fig. 4. Plot of L2/L'^ and Te/L'^, represented by Q ™d A respectively, versus 
l^-'i^h = L~^l'^ . Error bars are smaller than the size of the symbols. The dashed 
lines are simply to guide the eye. 



the honeycomb-lattice loop model defined by (11. 3p displays simultaneously the 
universal properties of a densely-packed loop model with loop fugacity n and 
those of a model with central charge c = 1 and thermal dimension = 1. The 
zero-temperature triangular-lattice antiferromagnetic Ising model has c = 1 
and Xt = 1, and when n = 1 the densely-packed loop model is in the perco- 
lation universality class. We may therefore expect the n = 1 FPL model to 
display some of the critical behavior of percolation. 



In Fig. m we plot the data for L2/L?' and Tg/L^ versus L and in Fig. [5] 
we plot G2/ L"^ versus L~^^h'''\ 



The data for n; and Ci were fitted to the ansatz 



0{L) = Oo + OiV + O2L 



S/2 



(3.6) 



with the exponents yi and y2 fixed to —2 and —4 respectively. We estimated 
Co = 0.028836(2) for m and and Oq = 0.02620(3) for In Fig. E] we plot rii 
and Ci versus 



Finally, we fit the data for the dimensionless ratio Qg to (13.61) with fixed 
exponents yi = —2 and y2 = —4, and with an additional correction term 
proportional to L^^l""-^. We estimate Co = 1.0248(4). 
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Fig. 5. Plot of versus L 
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^-5/24 _ Error bars are smaller than the size 



of the symbols. The dashed lines are simply to guide the eye. 
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Fig. 6. Plot of ni and Q, represented by □ and Q respectively, versus L~^. Error 
bars are smaller than the size of the symbols. The dashed lines are simply to guide 
the eye. 



3.3 Dynamic data 



For any observable O, we define its autocorrelation function 



po{t) := m)om - {o)\ 



21 



where (■) denotes expectation with respect to the stationary distribution. We 
then define the corresponding exponential autocorrelation time 



Texp,© := limsup — — , (3.7) 
t^ioo - log \po{t)\ 

and integrated autocorrelation time 

oo 

T-int,o ■■=7^ Poit) ■ (3.8) 

t=-oo 

Typically, all observables O (except those that, for symmetry reasons, are 
"orthogonal" to the slowest mode) have the same exponential autocorrelation 
time, so Texp,© = ^exp- However, they may have very different amplitudes of 
"overlap" with this slowest mode; in particular, they may have very different 
values of the integrated autocorrelation time, which controls the efficiency of 
Monte Carlo simulations 130]. 



The autocorrelation times typically diverge as a critical point is approached, 
most often like r ~ where ^ is the spatial correlation length and z is a dy- 
namic exponent. This phenomenon is referred to as critical slowing-down [SSEO]- 
More precisely, we define dynamic critical exponents ^exp and ^int.o by 



cxp 

7 



Texp ~ ^ 

On a finite lattice at criticality, C, can here be replaced by L. 



(3.9) 



During the simulations we measured the observables (except for Dq) only when 
the chain visited the Eulerian subspace, roughly every Te ~ L'^-'^^h iterations, 
or hits, of Algorithm [H However, it is natural when defining z^xp and Zmt.o via 
(13. 9p to measure time in units of sweeps of the lattice, i.e. L"^ hits. Since one 
sweep takes of order L'^^'^ visits to the Eulerian subspace, in units of "visits 
to the Eulerian subspace" we have r ~ L^+'^^h_ 

For each observable O = Afi, Vq, C2, G2 we computed poit) and Tint,© from 
our simulation data using the standard estimators discussed in pO]. By far 
the slowest of these observables is Afi. In Fig. [7| we plot puM/'TintM) ^"^^ 
observe that the decay is very close to being a pure exponential, suggesting 
^exp ~ -Zint,A/'r We fitted the Tint,A/'; data to the ansatz 

Tint = a + 6L^'-+2^\ (3.10) 

which produced the estimate 

zint,Ni = 0.515(8), 
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Fig. 7. Autocorrelation function pj\fi{t) versus t/T\a\.,Mi - clear that pj\fi{t) decays 
almost as a pure exponential, suggesting Zexp ~ -^intX; • 

suggesting 

^exp = 0.515(8). 

Assuming ^exp = ^intXn long-time decay of the autocorrelation function 
for any observable O should behave like po(t) ~ exp(— t/rint,A/'J- However, 
it was observed in [T7j that for the standard worm algorithm simulating the 
critical Ising model on the square and simple cubic lattices, some observables 
can have quite unusual short-time dynamics. Indeed it was found that Vq 
decorrelated in 0(1) hits and a detailed investigation of pvo{t) was presented. 
This phenomenon in which some observables decorrelate on time scales much 
less than L^^^p has been dubbed critical speeding-up [36|20|17] . We have not 
performed a detailed investigation of the behavior of pvoit) here, however we 
note that Ti^t,Vo ~ 0.5, independent of L, showing clearly that Vq certainly 
exhibits critical speeding-up under the dynamics of Algorithm [H For £2 and 
Q2 the short-time decay of p(t) appears to be intermediate between that of Afi 
and Vq. To illustrate, in Fig. Owe plot pc2{'t/'^intM)- appears that pc(t) has 
a short-time decay on a time scale strictly less than L^™p . Similar behavior is 
observed for pg^it)- 



4 Discussion 

We have formulated a worm algorithm that correctly simulates the FPL model 
on the honeycomb lattice when n = 1. Furthermore, we have rigorously proved 
its validity by showing that the corresponding Markov chain is irreducible 
and has uniform stationary distribution. Using standard duality relations this 
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Fig. 8. Autocorrelation function pc2{t) versus t/Ti^t,Afr 



algorithm can also be used to simulate the zero-temperature triangular-lattice 
antiferromagnetic Ising model. 

We have tested this worm algorithm numerically and estimate Zexp = 0.515(8), 
which suggests that it suffers from only mild critical slowing down. We observe 
that the dynamics of the algorithm exhibits the multi-time-scale behavior 
observed in [17j. It would be interesting to to examine the dynamic behavior 
of observables other than Afi in more detail, along the lines presented in jTTj . 
but this we leave to future work. We also obtained some interesting results 
regarding the static behavior of the FPL model, notably that the face-size 
moments appear to be governed by the magnetic dimension for percolation. 
This is consistent with the argument in [22] that the FPL model for general 
n displays simultaneously the universal properties of a densely-packed loop 
model and those of a model with central charge c = 1 and thermal dimension 



Finally, we note that one could in principle simulate (11. 3p with n > 1 by in- 
corporating appropriate connectivity checking into the Metropolis acceptance 
probabilities, or by combining an n = 1 worm algorithm with a "Chayes- 
Machta coloring" as described in [S7] . 
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A Topological constraints on fully-packed subgraphs of the honey- 
comb lattice 

The following lemmas describe some topological constraints on fully-packed 
spanning subgraphs of the honeycomb lattice with periodic boundary condi- 
tions. They are completely independent of any considerations regarding worm 
algorithms. We make essential use of Lemma lA.ll in the proof of Proposi- 
tion [221 

Lemma A.l. // {A, u,v) & TZ and u ^ v then du{A) = dy{A) . 

Lemma A. 2. Let {A,u,v) G TZ, and let Hi, Ui, Di denote, respectively, the 
number of vacant horizontal edges, the number of occupied up-pointing vertical 
edges, and the number of occupied down-pointing vertical edges, in row i. If 
row i contains no defects we have Hi = Ui = Di. 

Proof. Fix a configuration [A, u, v) G TZ, and a row i which contains no defects. 
Full-packing then implies that every vertex in this row has degree 2, so that 
for every vacant horizontal edge, one of its endpoints must be adjacent to an 
occupied up-pointing vertical edge and its other endpoint must be adjacent 
to an occupied down-pointing vertical edge, so Hi < Ui and Hi < Di. See 
Fig lA.li Conversely, if m is a vertex in row i which is adjacent to an occupied 
down-pointing vertical edge then precisely one of its horizontal edges must be 
vacant, so Di < Hi and therefore Di = Hi. Similarly, if f is a vertex in row i 
which is adjacent to an occupied up-pointing vertical edge then precisely one 
of its horizontal edges must be vacant, so Ui < Hi and therefore Ui = H^. □ 

Proof of Lemma \A.1[ Let us first note that there do indeed exist configu- 
rations with du{A) = dv{A) when {A,u,v) G 71 and u ~ v. Indeed, if 
{A, u,u) E 71 and t> ~ m then {AAuv, u, v) E 7Z; if uv E A then du{AAuv) = 
1 = dy{AAuv), whereas ii uv ^ A then du{AAuv) = 3 = dy{AAuv). 

Suppose on the contrary that {A,u,v) G 7Z with -u ~ u, but du{A) ^ dy{A). 
Since u v we must have du{A),dy{A) G {1,3}, so that either du{A) = 3 
and dy{A) = 1, or vice versa. Let us assume (without loss of generality) the 
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Fig. A.l. If the horizontal edge uv is vacant, and neither u nor w is a defect, the 
remaining edges incident to both u and v are forced to be occupied, in a fully 
packed configuration. Conversely, if the up-pointing vertical edge vw is occupied, 
then precisely one of the horizontal edges incident to v must be vacant (here uv), 
and the down-pointing vertical edge ut must then be occupied. In the diagram, thick 
edges are occupied, thin edges are vacant, and dotted edges are unconstrained by 
the state of uv. 




Fig. A. 2. Neighborhood of the horizontal edge uv with da{A) = 3 and dv{A) = 1. 




former. There are two possibilities for the edge uv; either uv is a horizontal 
edge, so that u and v lie in the same row, or uv is a vertical edge, so that u 
and V lie in adjacent rows. 

Suppose uv is a horizontal edge lying in row i, denote f 's other horizontal 
edge by vw, and suppose that the vertical edge uu' is up-pointing, so that the 
vertical edge vv' must be down-pointing. See Fig. IA.2[ The up-pointing vertical 
edges uu' and ww' are both occupied. Suppose there are n other occupied up- 
pointing vertical edges incident to row i, so there are n + 2 in total. Each 
of these other n occupied up-pointing vertical edges is incident to a degree 2 
vertex in row i. Let a be such a vertex, then a must have one of its horizontal 
edges vacant, call it ab. By assumption we have a ^ u,w, so that b ^ v, 
and so db{A) = 2 and b must have its vertical edge (which is down-pointing) 
occupied. Therefore, every one of the n occupied up-pointing vertical edges 
other than uu' and ww' corresponds to an occupied down-pointing vertical 
edge. Conversely, if there is an occupied down-pointing vertical edge incident 
to some b ^ V in row i then b must have one of its horizontal edges vacant, 
call it ab. Since b ^ v and ab is vacant we have a ^ u,w, so that da{A) = 2. 
Therefore a must have its vertical edge (which is up-pointing) occupied, and 
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this edge is neither uu' nor ww' . Therefore there are n occupied down-pointing 
and n + 2 occupied up-pointing vertical edges incident to row i. Now, since 
no other row contains a defect, Lemma [A. 21 tells us that all rows below row 
i will have n occupied up-pointing and down-pointing vertical edges, whereas 
all rows above row i will have n + 2 occupied up-pointing and down-pointing 
vertical edges. However, it is impossible for this to occur if we have periodic 
boundary conditions, and so we have a contradiction. Of course, if we assume 
instead that u is down-pointing and v up-pointing then an entirely similar 
argument leads to a similar contradiction. Therefore if {A, u,v) G 7^ and uv 
is a horizontal edge we must have du{A) = dy{A). 

The converse situation where uv is a vertical edge can be treated in a similar 
manner. We omit the details. □ 



B Lemmas used in the proof of Proposition 12.31 

Lemma B.l. Let {A,v,v) G S with vv' ^ A. Then whenever dy{A) = 2 we 
have 

(A, V, v) ^ {AU vv', v', v) (A, v', v'), 
and whenever d^^A) = dyi{A) = 2 we have 

{A,v,v)^{A\vv",v",v)-^{A,v",v"), 

for both v" ~ V with v" ^ v' . 

Lemma B.2. Let H G denote the fully-packed configuration in which 
every horizontal edge is occupied and every vertical edge is vacant. Suppose 
{A,u,u) G S with du{A) = 2 and A ^ H. Then there always exists a vacant 
horizontal edge vv' for which {A, u, u) — > {A, v, v). 

Lemma B.3. If {A, v,v) eTZ then {A, v, v) ^ {A, u, u) for any pair of vertices 
u and V. 

Lemma B.4. // {A, v,v) eTZ and -u ~ f then 

(y4,t>,t>) <-> {A/\uv,u,v). 
Lemma B.5. Let {A,u,v) G TZ with du{A) = 1 and suppose uu' ^ A. Then 

{A, u, f ) (A U uu', u' , v) 
Lemma B.6. // [A, u,v) with du{A) = 3 then for each u' ^ u 

{A, u, v) ^ {A \ uu', u', v). 
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Proof of Lemma \B. 1[ For any {A, v,v) E S with vv' ^ A we have 
P:^[{A,v,v)-^{AUvv',v',v)] = 1/6. 

Furthermore, if dy{A) = 2 then d^{A U vv') = 3 so 

P:^[iAUvv',v',v)^{A,v',v')] = l/6 

and we have 

{A, V, f ) — > (A U vv', v', v) — > {A, v', v'). 

If in fact dv{A) = 2 = dvi{A) then for both v" ~ v with v" ^ v' we have 

P'^[{A,v,v)^{AVJvv',v,v')] = l/Q 
P'Jy{A U vv', V, v') (A U W \ vv", v", v')] = 1/6 
[(A U w' \ vv", v", v') ^{A\ vv", v", v)] = 1/6 
P:^[{A\vv",v",v)^{A,v",v")] = 1/4 



so that 



(A, V, v)-^iA\ vv", v", v) (A, v", v"). 



□ 



Proof of Lemma \B.^ Denote by uu' ^ A the unique vacant edge incident to 
u. There are two possibihties: either du'{A) = or du'{A) = 2. If du'{A) = 
then u' has both its incident horizontal edges vacant, and since Lemma IB.ll 
implies {A, u, u) — >■ {A, u' , u') there is nothing more to show. If on the other 
hand dui{A) = 2 then Lemma fB . 1 1 implies {A, u, u) {A, v, v) for every v ^ u. 
If any of the {A, v, v) have a vacant horizontal edge incident to v we are done. 
Otherwise we re-apply Lemma [B. II to {A, v, v) for every v ~ u. In this way we 
must eventually arrive at some {A, w, w) for which there is a vacant horizontal 
edge incident to w. Transitivity implies {A, u, u) — > {A, w, w) and the stated 
result follows. □ 



Proof of Lemma \B. 3[ If {A,v,v) G TZ then in fact {A,u,u) G TZ for any u. 
Since every vertex has degree 2 we can apply Lemma [B.ll to {A,v,v) to see 
that {A, V, v) {A, v', v') for any v' ~ v, but we can equally apply it to 
{A,v',v') to see that {A,v',v') {A,v,v). So for any pair of neighboring 
vertices w ~ f ' we have {A,v,v) {A,v',v'). Since the lattice is connected 
and every vertex has degree 2 this immediately extends, via transitivity of 
to [A, V, v) ^ {A, u, u) for any arbitrary pair or vertices u, v. □ 



Proof of Lemma B.4 Let {A,v,v) G TZ and m ~ f . Lemma [B.ll immediately 
implies 

{A,v,v) {AAuv,u,v) 



28 



and combining Lemma IB.ll with Lemma IB. 31 we further obtain 

{AAuv,u,v) — > {A,u,u) {A,v,v). 
Therefore {A,v,v) {AAuv,u,v). □ 

Proof of Lemma \B.^ Suppose that u' ^ v. Then 

P'^ [{A, u, v)^{A\^ W, u', t;)] = i 
and since u' ^ v imphes du'^A) = 2 we have du'^A U uu') = 3, so that 

P^[(A U uu', u\ v) (A u, v)] = ^. (B.l) 
Therefore (A, m, t>) ^ (A U uu' , u', v) when m' 7^ w. 

Conversely, suppose u' = v. Then Lemma [B.4I imphes that {AUuv,v,v) ^ 
{A,u,v). □ 

Proof of Lemma \B.6[ Suppose that u' 7^ v. Since = 3 

P:^[{A,u,v) ^ {A\uu',u',v)] = ^. 
Furthermore, since u' ^ v we have du>{A) = 2 and hence du>{A \ uu') = 1, so 

P:^[{A\uu',u',v)^{A,u,v)] = ^. 
Therefore {A, u, v) ^ {A \ uu', u', v) when u' 7^ v. 

Conversely, suppose u' = v. Then Lemma lB.41 implies that [A \ uv,v,v) ^ 
{A,u,v). □ 

References 

[1] G. H. Wannier, Phys. Rev. 79 (1950) 357. 

[2] J. Stephenson, J. Math. Phys. 5 (1964) 1009. 

[3] H. J. H. B. Nienhuis, H. W. J. Blote, J. Phys. A 17 (1984) 3559. 

[4] H. W. J. Blote, M. P. Nightingale, Phys. Rev. B 47 (1993) 15046. 

[5] S. L. A. de Quieroz, E. Domany Phys. Rev. E 52 (1994) 4768. 

[6] X. F. Qian, H. W. J. Blote, Phys. Rev. E 70 (2004) 036112. 



29 



[7] X. F. Qian, M. Wegewijs, H. W. J. Blote, Phys. Rev. E 69 (2004) 036127. 

[8] E. Rastclli, S. Rcgina, A. Tassi, Phys. Rev. B 71 (2005) 174406. 

[9] K. Xia, X.-Y. Yao, J.-M. Liu, Frontier of Phys. in China 2 (2007) 191. 

[10] G. M. Zhang, C. Z. Yang, Phys. Rev. B 50 (1994) 12546-12549. 

[11] P. D. Coddington, L. Han, Phys. Rev. B 50 (1994) 3058. 

[12] D. Kandel, R. Ben-Av, E. Domany, Phys. Rev. B 45 (1992) 4700. 

[13] A. Dhar, P.Chaudhuri, C. Dasgupta, Phys. Rev. B 61 (2000) 6227. 

[14] R. H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 58 (1987) 86-88. 

[15] N. Prokof'ev, B. Svistunov, Phys. Rev. Lett. 87 (2001) 160601. 

[16] M. Jerrum, A. Sinclair, SIAM J. Comput. 22 (1993) 1087. 

[17] Y. Deng, T. M. Garoni, A. D. Sokal, Phys. Rev. Lett. 99 (2007) 110601. 

[18] U. Wolff, Simulating the All-Order Strong Coupling Expansion L Ising Model 
Demo, arXiv:0808.3934vl [hep-lat] . 

[19] M. Sweeny, Phys. Rev. B 27 (1983) 4445-4455. 

[20] Y. Deng, T. M. Garoni, A. D. Sokal, Phys. Rev. Lett. 98 (2007) 230602. 
[21] L. Chayes, J. Machta, Physica A. 254 (1998) 477-516. 

[22] B. Nienhuis, Phys. Rev. Lett. 49 (1982) 1062. 
[23] B. Nienhuis, J. Stat. Phys. 34 (1984) 731. 

[24] P. Di Francesco, P. Mathieu, D. Senechal, Conformal Field Theory, Springer- 
Verlag, New York, 1997. 

[25] O. Schramm, Israel J. Math. 118 (2000) 221. 

[26] S. Rohde, O. Schramm, Ann. Math. 161 (2005) 883. 

[27] G. Lawler, Conformally Invariant Processes in the Plane, American 
Mathematical Society, Providence, 2005. 

[28] H. Blote, B. Nienhuis, Physica A. 160 (1989) 121. 

[29] H. W. J. B16te, B. Nienhuis, Phys. Rev. Lett. 72 (1994) 1372. 

[30] A. D. Sokal, Monte Carlo methods in statistical mechanics: Foundations and 
new algorithms, in: P. C. C. DeWitt-Morette, A. Folacci (Eds.), Functional 
Integration: Basics and Applications, Plenum, New York, 1997, pp. 131-192. 

[31] J. G. Kemeny, J. L. Snell, Finite Markov Chains, Springer- Verlag, New York, 
1976. 

[32] J.-S. Wang, Phys. Rev. E 72 (2005) 036706. 



30 



[33] M. losifescu, Finite Markov Processes and Their Applications, John Wiley & 
Sons, Bucharest, 1980. 

[34] G. Grimmett, D. Stirzaker, Probability and Random Processes, 3rd Edition, 
Oxford University Press, New York, 2006. 

[35] P. C. Hohenberg, B. I. Halperin, Rev. Mod. Phys. 49 (1977) 435. 

[36] Y. Deng, T. M. Garoni, A. D. Sokal, Phys. Rev. Lett. 98 (2007) 030602. 

[37] Y. Deng, T. M. Garoni, W. Guo, H. W. J. Blote, A. D. Sokal, Phys. Rev. Lett. 
98 (2007) 120601. 



31 



